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I. INTRODUCTION 

Within cosmology, astrophysics or geophysics one often 
needs to deal with electrically conducting plasmas at high 
kinematic and magnetic Reynolds numbers where mag- 
netic fields are dynamically important. Indeed, much of 
the turbulence in the interstellar medium is magnetohy- 
drodynamic in nature. 

Hydromagnetic turbulence has been explored exten- 
sively in connection with the generation of large scale 
magnetic fields in astrophysical bodies such as planets, 
stars, accretion discs and galaxies through dynamo the- 
ories. Non-driven, freely decaying turbulence may also 
be of interest in connection with both the physics of the 
interstellar medium and cosmology. Our interest was in- 
spired by the cosmology of primordial magnetic fields, 
which is sometimes considered as a possible source for 
providing the seed for the galactic dynamo [fij . 

There have been various related works on decaying 
MHD turbulence, by authors interested in different con- 
texts H -fjj. Most directly comparable to our work, 
Biskamp and Miiller Q studied the energy decay in 
incompressible 3D magnetohydrodynamic turbulence in 
numerical simulations at relatively high Reynolds num- 
ber, and in a companion letter Q studied the scaling 
properties of the energy power spectrum. 

We are here especially interested in the inverse cascade 
of magnetic helicity whereby magnetic energy is trans- 
ferred from small to large scale fluctuations. This is im- 
portant for a primordial magnetic field to reach a large 
enough scale with sufficient amplitude to be relevant for 
seeding the galactic dynamo ||] . 

It should be noted that due to the conformal invari- 
ance of MHD in the radiation era the MHD equations 
in an expanding universe can be converted into the rela- 
tivistic MHD equations in flat spacetime by an appropri- 
ate scaling of the variables and by using conformal time 
H . The equations of (9) differ slightly from the ordinary 
non-relativistic MHD equations. However, in order to 
facilitate comparison with earlier work, we use the non- 



relativistic equations. 

We perform 3D simulations both with and without 
magnetic helicity, starting from statistically homoge- 
neous and isotropic random initial conditions, with power 
spectra suggested by cosmological applications. We find 
a strong inverse cascade in the helical case, with equivo- 
cal evidence for a weak inverse cascade when only helicity 
fluctuations are present. In the helical case we also find a 
self-similar power spectrum with an approximately /c~ 2,5 
behavior at high k. We present energy decay laws which 
are comparable to those found in the incompressible case 
by Biskamp and Miiller ||, and in the compressible case 
by Mac Low et al. 



II. THE MODEL 

We consider the equations for an isothermal compress- 
ible gas with a magnetic field, which is governed by the 
momentum equation, the continuity equation, and the 
induction equation, written here in the form 



du 

~dt 



u ■ Vu - c 2 Vlnp- 



J x B 



<91np 
dt 



.t fv 2 u+iw-u 



= u • V In p — V • u, 



dA 

— = uxB + 77V 2 A, 
at 



(1) 



(2) 



(3) 



where B = V x A is the magnetic field in terms of the 
magnetic vector potential A, u is the velocity, J is the 
current density, p is the density, p is the dynamical vis- 
cosity, and r\ is the magnetic diffusivity. 

The code for solving these equations |T(]] uses a vari- 
able third order Runge-Kutta timestep and sixth order 
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explicit centered derivatives in space. All our runs are 
performed on a 120 3 grid, and we use periodic boundary 
conditions, which means that the average plasma density 
(po) = Po is conserved during runs. Here po is the value 
of the initially uniform density, and the brackets denote 
volume average. 

We adopt nondimensional quantities by measuring u 
in units of c, where c is the speed of light, k in units of k\, 
where k\ is the smallest wavenumber in the box, which 
has a size of Lbox = 27r, density in units of po = 1, and B 
is measured in units of y/p p c, where po is the magnetic 
permeability. This is equivalent to putting c = fci = 
Pa = Pa = 1- In the following we will refer to the mean 
kinematic viscosity v which we define as v = p/po- The 
sound speed c s takes the value c s = l/y/3, as appropriate 
for a relativistic fluid. With c = 1, the unit of time is 
such that the light crossing time of the box is 2w. 

Our equations are similar to those for the relativistic 
gas in the early universe using scaled variables and con- 
formal time for non-relativistic bulk velocities ||. We 
expect our results to change little using the true rela- 
tivistic equations, as our advection velocity is at most 
only mildly relativistic, and this only at the beginning of 
the simulation. 



III. ON THE ROLE OF THE INVERSE CASCADE 



The magnetic helicity Hm is given by 



H M = I A • B 



(4) 



and characterizes the linkage between magnetic field 
lines. Hm is conserved in the absence of ohmic dissipa- 
tion, although it is still possible to have local, small scale 
helicity fluctuations. Helicity plays an important role in 
dynamo theory where turbulence is driven. 

In many astrophysical and cosmological situations the 
magnetic Reynolds number ReM is very large. We define 
the magnetic Reynolds number as ReM = Lv/rj, where 
L and v are the typical length scale and velocity of the 
system under consideration and 77 is the resistivity. The 
magnetic Reynolds number is a measure of the relative 
importance of flux freezing versus resistive diffusion. In a 
cosmological context this number can be extraordinarily 
large: causality imposes the weak limit L < ct and rel- 
ativity demands v < c. With conductivities relevant to 
the era when the electroweak phase transition took place 
fl3| ], one can in principle obtain a magnetic Reynolds 
number of about 10 16 . This is often taken to mean that 
the magnetic field is frozen into the plasma, and the scale 
length of the field increases only with the expansion of 
the Universe. 

However, this simple picture does not necessarily give a 
full description of the dynamics because the MHD equa- 
tions, especially at high Reynolds numbers where non- 
linear terms are important, exhibit turbulent behavior, 



which can lead to a redistribution of magnetic energy 
over different length scales ||. Energy in a turbulent 
magnetic field can undergo an inverse cascade and be 
transferred from high frequency modes to low frequency 
modes, increasing the overall comoving correlation length 
|]lT|| . This process is due to the nonlinear terms giving 
rise to interactions between many different length scales. 

We will take the initial primordial power spectrum as 
given and address the question of how such a primor- 
dial spectrum evolves as a consequence of the nonlinear 
equations of motion. 



IV. INITIAL CONDITIONS 

Since one of the aims of the present work is to inves- 
tigate the role of magnetic helicity in the inverse cas- 
cade we describe how the initial conditions for our sim- 
ulations were set up. We chose our initial condition by 
setting up magnetic fluctuations with an initial power 
spectrum PM{k) = (B£ • Bk) ~ fc" in Fourier space (and 
averaged over shells of constant fc = |k|), for low val- 
ues of the wavenumber fc, using a exponential cutoff k c . 
[The shell-averaged power spectrum, PM{k), is not to 
be confused with the shell-integrated energy spectrum, 
Em = 47rfc 2 x ^PM(k), which is shown in the plots be- 
low.] 

The magnetic field fluctuations are drawn from a Gaus- 
sian random field distribution fully determined by its 
power spectrum in Fourier space according to the follow- 
ing procedure. For each grid point we use the correspond- 
ing wavenumber to select an amplitude from a Gaussian 
distribution centered on zero and with the width 



Pm(k) =P M ,ofc n exp(-(fc/fc c ) 4 ) 



(5) 



where k = |k|. We then transform the field back into real 
space to obtain the field at each grid point. This is done 
independently for each field component. 

There is a requirement in cosmology that n > 2, which 
is set by causality demanding that the correlation func- 
tion of the magnetic field vanishes at large distances, and 
the fact that the magnetic field is divergence- free [fT4| . 
In the simulations presented we chose the slope of the 
power spectrum to be n — 2. We also chose k c = 30, 
unless specified otherwise, which gives a power spectrum 
peaked at a relatively large value of k. Biskamp and 
Miillcr j(||7j started with a spectrum peaked at k c = 4, 
which may account for the different slope in the late- time 
power spectrum which we observe (see Section V A ). 

Our velocity power spectrum was chosen in a similar 
way, but with n = corresponding to white noise at large 
scales (there is no requirement for incompressibility in the 
early Universe). The initial magnetic energy was taken 
equal to the kinetic energy, and had the value 5 x 1CP 3 
in all runs, as the primordial field is thought likely to be 
weak. 



2 



In order to introduce a non-zero average magnetic he- 
licity into the system it is useful to represent the vector 
potential in terms of its projection onto an orthogonal 
basis formed by e+, e and k. The two basis vectors e + 
and e_ can be chosen to be the unit vectors for circular 
polarization, right-handed and left-handed respectively. 
That is e± = ei ± ie.i where ei and e.i are unit vectors 
orthogonal to each other and to k. They are given by 
ei = k x z/|k x (k x z)| and e 2 = k x (k x z)/|k x z| 
respectively, z is a reference direction. 

Note that since 

ik x e s = ske s (6) 

where s = ±1, this corresponds to an expansion of the 
magnetic vector potential into helical modes. 

Using these basis vectors it is easily seen that the mag- 
netic energy spectrum is 

E M (k) = 2^fc 2 (|B k | 2 ) (7) 

where the amplitude of the magnetic field is given by 

|B k | 2 = (|A+| 2 + |A k | 2 )|k| 2 (8) 

and the expression for the magnetic helicity spectrum 
H M (k) is 

H M (k) = 4^fc 2 (A* B k ) (9) 

where 

A k -B k = (L4+| 2 -L4 k | 2 )|k|. (10) 

The function Hyi{k) is a sensitive measure of the cor- 
relation between the vector potential and the magnetic 
field. Hyi{k) may, of course, be positive in one part of 
Fourier space and negative in another part. It is, how- 
ever, bounded in magnitude by the inequality 

\H u (k)\ <2k- 1 E M (k). (11) 

A field which saturates the above inequality is maximally 
helical. 

The amplitudes can be chosen independently, pro- 
vided A*J^ = A k , which is just the condition that the 
vector potential be real. Therefore it is possible to adjust 
the amplitudes \A£\ and \A^\ freely and in so doing ob- 
taining a magnetic field with arbitrary magnetic helicity. 
With our method we are able to put statistically random 
but maximally helical fields in our initial conditions. In 
our runs with initial helicity we take Hm = H max . 

Because we evolve our dynamical fields on a discrete 
lattice we have to be careful when using derivative opera- 
tions in Fourier space. In general, the wave vector, which 
is an eigenvalue of the derivative operator, needs to be 
replaced by some function k e s(k), which is an eigenvalue 
of the discrete derivative operator on the lattice. In our 
case we have, for the sixth order explicit centered deriva- 
tive 



k eS (k) = ± [sin(3fe) - 9 sin(2fc) + 45 sin(fc)] . (12) 

In order to be consistent with the scheme used in the 
simulation, we use k e {f(k) when calculating the initial 
condition in Fourier space. 

V. RESULTS 

In all runs the mean kinematic viscosity v and the re- 
sistivity r\ were chosen to be equal with values between 
v = t) = 5 x 10~ 4 — 5 x 10~ 5 . In our simulations we typi- 
cally obtain Reynolds numbers of the order of 100 — 200. 
The Reynolds numbers in our simulations are evaluated 
using the magnetic Taylor microscale which we calcu- 
late here as the ratio of the rms magnetic field and the 
rms current density, Lt = 2TtB ims / J rros . The 2ir factor 
is here included so that Lt represents the typical wave 
length (and not the inverse wavenumber) of structures in 
the current density. 

A. Spectral evolution 

The inverse magnetic cascade for decaying MHD tur- 
bulence is best visualized in terms of magnetic energy 
spectra Ey^(k) because information on nonlinear interac- 
tion between different scales is contained in Eyi(k). In 
Fig. [l] we show a run with initial magnetic helicity. 
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FIG. 1. Magnetic energy spectrum Eyi(k) for a run with 
finite magnetic helicity. v — rj — 5 x 1CF J . The times shown 
are 0, 1.0, 4.6, 10.0, 21.5 and 46.3. The initial spectrum is in- 
dicated by the dashed line. At low wavenumbers k the energy 
spectrum £m(&) increases with time. 

In Fig. ^] we see evidence for a dual energy transfer both 
toward higher and lower wavenumbers. The inverse cas- 
cade is characterized by the transfer of energy from small 
scale structures in the magnetic field to larger ones. In 
Fig. [I] this behavior is clearly seen as indicated by the 
rise in the energy spectrum at small wavenumbers. Some 
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energy is also being transported to smaller scales where 
the spectrum is decaying due to diffusive effects. We also 
note that at wavenumbers above the peak k p (t) the spec- 
trum develops a power law shape. This power law has 
approximately a fc -2 5 slope. This differs from the ap- 
proximately fc -5 / 3 law found by Miiller and Biskamp [Q. 
We suggest that this is due to finite size effects, which af- 
fect the spectrum if the initial scale separation between 
k p and the smallest wavenumber in the box (k = 1) is 
insufficient, and if the flow is strongly helical so that its 
spectrum is governed by inverse cascading. In order to 
check this we have performed a run with larger initial 
length scale, k c = 5. In this case the magnetic energy 
spectrum develop into an approximate k~ 5 ^ 3 law at late 
times. However, this occurs only after the peak of the 
spectrum has left the simulation box, i.e. after finite size 
effects have begun to play a role. 

To check if the magnetic field evolution is self-similar 
one can make the following ansatz for the energy spec- 
trum 

E M (k,t) = t(ty q m(kO. (i3) 

Here £ is the characteristic length scale of the magnetic 
field, taken to be the magnetic Taylor microscale defined 
above, and q is a parameter whose value is some real 
number. We call gu the magnetic scaling function. 
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FIG. 2. The ma gne tic scaling function givi(fc£) described in 
the text, equation Q13) , versus The straight lines indicate 
the power laws tx (fe£) 4 ' an( l K (&£)~ 2 5 respectively. 

In Fig. ||we have plotted £(t) q EM(k,t) versus the scaled 
variable fc£ (t) . The value of the parameter q in this run is 
q = 0.7. It is seen that for each different value of time t, 
the data collapses onto a single curve given by the scaling 
function <7m demonstrating the self-similarity of the 
magnetic field evolution. 

We also performed runs in which the magnetic helicity 
was zero, in the statistical sense. Magnetic helicity was 
present due to fluctuations, but was of very small am- 
plitude. In these runs no significant inverse cascade was 
observed. Fig. || shows the energy spectrum for such a 



run with only small magnetic helicity fluctuations present 
in the initial conditions. It is seen that only a weak in- 
verse cascade is present at the lowest wavenumbers, much 
smaller than in the helical case. However, that it seem 
to be present at all is interesting as the effect could be- 
come more pronounced for higher Reynolds numbers. It 
is possible that this effect is due to the magnetic helicity 
fluctuations even though they were small. One simula- 
tion was performed with identically zero initial magnetic 
helicity fluctuations. In this case random fluctuations de- 
velop rapidly and no differences between the two cases, 
were observed. 
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FIG. 3. Magnetic energy spectrum Eu(k) for a run with no 
net magnetic helicity. v = r) = 1 x 10~ 4 . Here k c = 10. The 
times shown are 0,2.2,4.6,10.0,21.5 and 46.3. The initial 
spectrum is indicated by the dashed line. The peak of the 
energy spectrum _Em(&) is decreasing with increasing time. 

B. Energy decay 

In Fig. | we show the time evolution of the magnetic 
energy EM{t) and the kinetic energy E^{t) for a run with 
initial helicity and a fc 4 initial energy spectrum slope. 
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FIG. 4. Time evolution of the magnetic energy Eyi(t) and 
the kinetic energy Escif) in the case where there is initial 
magnetic helicity. u — r\ — 5 x 1CP 5 . The straight lines 
indicate the power laws oc t~ ' 7 and oc t~ respectively. 

It is seen that the asymptotic decay rate for -Em(^) is 
approximately £~ 0,7 . The Reynolds number for this run 
was around Re ~ 200 at late times. In another run with 
Re ~ 100 the decay rate was seen to be t~°' s , so there 
seems to be a dependence of the decay rate of the mag- 
netic field on the Reynolds number and perhaps the re- 
sulting slope of the spectrum. 

The kinetic energy also decays with a power law behav- 
ior at late times. In the case of runs with initial helicity 
the kinetic energy -Ek(^) decays with a different, faster 
rate than Eu(t)- The asymptotic decay rate is close to 
t . In runs without initial helicity the decay rates of 
Eyi(t) and E^(t) are approximately the same, close to 
f~ lA . 

In our runs with Ek = Em initially, the kinetic en- 
ergy spectrum shows no evidence of an inverse cascade at 
any scale. However, when the initial velocity distribution 
is zero the kinetic spectrum grows on all scales initially 
and in the low wavenumber region the energy continues 
to grow even after the high wavenumber modes start to 
decay. 



C. Coherence length evolution 

During the course of the simulations the initially small 
scale structures gain in size. A convenient length scale 
is the magnetic Taylor microscale Lt, which was defined 
above. This length scale is mostly characteristic of the 
small scales but even they grow during the course of the 
simulations. 

In Fig. [s] we show the evolution of for a run with 
initial helicity. 
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FIG. 5. Time evolution of the magnetic Taylor microscale 
for the case with initial magnetic helicity. v = tj — 5 x 1CP 5 . 
The straight line indicates the power law oc t ' 5 . 



The asymptotic behavior of the length scale is seen to 
grow approximately as Lt ~ i ' 5 - 

In runs with non-helical initial conditions the growth 
of the magnetic Taylor microscale is slower than in the 
case of helical initial conditions. In this case the magnetic 
Taylor microscale grows approximately as Lt ~ t 0A . 

The discussion so far has mainly been concerned with 
the evolution of causally generated magnetic fields using 
an initial fc 4 slope in the magnetic energy spectrum. Now 
we briefly comment on the other cases we have looked at. 
For a white noise initial spectrum Eyi{k) ~ k 2 , the evo- 
lution is qualitatively and quantitatively similar to the 
causal case. For helical fields we observe an inverse cas- 
cade, while for non-helical fields a much smaller inverse 
cascade is present only for the lowest modes. 



VI. DISCUSSION 

Our simulations show the decay rate of magnetic en- 
ergy for compressible turbulence being sensitive to the 
initial helicity of the magnetic field configuration. A sim- 
ilar result was found in in the case of incompressible 
turbulence. The fact that magnetic helicity is conserved 
(except for resistive changes), and the magnetic energy 
decays slower for helical fields, is connected with the ob- 
served inverse cascade in which magnetic energy is trans- 
ported toward larger scales because of nonlinear dynam- 
ics. 

The decay of kinetic energy does not seem to depend 
on the initial helicity and its decay rate Exit) ~ i -1 ' 1 
is consistent with the earlier work of H§|. Note that in 
the helical case we observe the kinetic energy decaying 
more rapidly than the magnetic one; this behavior was 
also found in ||. 

While these results are not directly applicable to the 
evolution of primordial magnetic fields in the early uni- 
verse, they do suggest that nonlinear magnetohydrody- 
namical effects may play an important role in this case. 

In any case, it is interesting to compare our results 
with the work of other authors interested in the decay 
properties of cosmological mag netic fields Ideal 
MHD has a scale invariance which leads to the scaling 
law H0 



E M (t,k) 



-l-2h 



(14) 



where ip is an unknown function, related to gyi- Assum- 
ing it is peaked somewhere, and h < 0, the characteristic 
scale of the field goes as L(t) ~ It is also of- 

ten assumed that ^(0) exists and is non-zero: thus h is 
determined by the initial power spectrum. Hence for a 
magnetic power spectrum of index n, h = — (n+3)/2 and 



(15) 



This law can also be recovered by assuming that the char- 
acteristic time scale for the decay of turbulence on a scale 
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I is the eddy turnover time tj = l/vu where vi ~ ;~( n + 3 ) 
is the velocity averaged on a scale I |16| . If the character- 
istic scale of the field is that scale which is just decaying, 
then tl ~ t, and we again find Eq. (|l5|). One should 
note that these arguments ignore helicity conservation. 

We recall that our non-helical runs had n — 2 for the 
magnetic power spectrum and n = for the velocity 
power spectrum. The observed growth law for the mag- 
netic Taylor microscale, t 0A , is not consistent with the 
predicted power law for n = 2, although it does square 
with the growth law for n — 0, and it is possible that 
the growth in the magnetic field length scale is being 
controlled by the velocity field. Simulations at higher 
Reynolds numbers seem required to resolve this issue. 

One would expect on integrating the helicity power 
spectrum that Hm ~ LjEm, where L% is the integral 
scale. We would expect that Lj ~ Lt, and hence, if 
magnetic helicity is conserved, 

E M ~ L^}. (16) 

However, magnetic helicity is not conserved exactly: we 
observe a decrease in Hm by a factor of about 2 in a run 
with viscosity v = 5 x 10~ 5 . Indeed, with Lt ~ t - 5 wc 
find a somewhat steeper relation: Em ~ t~ - 7 ~ L^, 1A . 

Finally, it is interesting to note that Son's numerical 
simulations of decaying turbulence [ pl| , performed in the 
Eddy-Damped Quasi-Normal Markovian (EDQNM) ap- 
proximation, show some evidence of a power law devel- 
oping at high k, the slope being close to fc -2 5 , although 
there was no net helicity present, and no inverse cas- 
cade. Furthermore, Field and Carroll |T^| , again using 
the EDQNM approximation, found that there were self- 
similar solutions with Em ~ t~ 2 / 3 ~ L^, 1 . 

VII. CONCLUSIONS 

We have shown that for an isothermal and compress- 
ible magnetized turbulent fluid, when undergoing a pro- 
cess of free decay, a substantial inverse cascade is present 
for helical magnetic field configurations, which transfer 
energy from smaller scale magnetic fluctuation to larger 
scale ones. For non-helical magnetic fields only a weak 
inverse cascade was observed on the largest scales. 

The energy spectrum of the magnetic field shows evi- 
dence for a self-similar evolution with a development of 
a power law of roughly fc -2 5 beyond the peak. Decay 
laws for both the kinematic and magnetic energy were 
found. The kinetic energy decay was approximately t~ 1A 
for both helical and non-helical magnetic fields. The de- 
cay of the magnetic field energy was found to be strongly 
dependent on the the initial helicity, decaying roughly 
as t~ 0J and t -1-1 for helical and non-helical initial con- 
ditions respectively. For the helical case, the magnetic 
energy decay rate showed a dependence on the Reynolds 
number, with a slower decay rate for larger Reynolds 
numbers. 



We also observed power law behavior in the charac- 
teristic length scale of the magnetic field, defined as the 
Taylor microscale Lt- In the helical case Lt ~ i ' 5 , 
whereas for non-helical fields the growth was somewhat 
slower, Lt ~ t 0A , and we ascribe the faster growth rate 
to the presence of the inverse cascade in the helical case. 
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